# Импротируем необходимые функции
import numpy as np # работа с массивами и линейной алгеброй
import matplotlib.pyplot as plt # для отрисовки графиков
import pandas as pd # для чтения и работы с данными
from scipy.optimize import curve_fit # фитирующая процедура
def psi_f(x):
return np.exp(-x**2/2)/np.sqrt(np.sqrt(np.pi))
def integrate_hist(hist, bin_edges):
integral = 0
for i in range(len(hist)):
integral += hist[i]*(bin_edges[i+1] - bin_edges[i])
return integral
def distr(df):
pedestrians, bin_edges = np.histogram(df, bins = 1000)
#print bin_edges
width = bin_edges[1] - bin_edges[0]
plt.figure(figsize=(7,7))
plt.bar(bin_edges[:-1], pedestrians, width = width)
plt.xlabel('Ось x'.decode('utf-8'))
plt.ylabel('Количество пешеходов (N)'.decode('utf-8'))
plt.title('Распределение пешеходов по оси'.decode('utf-8'))
plt.legend()
plt.show()
pedestrians_sq = pedestrians*pedestrians
norm = integrate_hist(pedestrians_sq, bin_edges)
psi_sq = pedestrians_sq/norm
psi = pedestrians/np.sqrt(norm)
plt.figure(figsize=(7,7))
plt.bar(bin_edges[:-1], psi, width = width)
plt.xlabel('x'.decode('utf-8'))
plt.ylabel('$\psi(x)$'.decode('utf-8'))
plt.title(''.decode('utf-8'))
###
x = np.arange(-4, 4, 0.1)
plt.plot(x, psi_f(x), color = 'black', label = '$\psi(x)$, теоретическая'.decode('utf-8'))
plt.legend()
plt.show()
plt.figure(figsize=(7,7))
plt.bar(bin_edges[:-1], psi_sq, width = width)
plt.xlabel('x'.decode('utf-8'))
plt.ylabel('$\psi^2(x)$'.decode('utf-8'))
plt.title('Плотность вероятности'.decode('utf-8'))
plt.plot(x, psi_f(x)**2, color = 'black', label = '$\psi^2(x)$, теоретическая'.decode('utf-8'))
plt.legend()
plt.show()
def drawNE(df):
plt.figure(figsize=(20,7))
plt.scatter(np.arange(df.shape[0]), df['N'], linewidths=0.1, label = 'Количество пешеходов'.decode('utf-8') )
plt.xlabel('Номер итерации'.decode('utf-8'))
plt.ylabel('Количество пешеходов (N)'.decode('utf-8'))
plt.title('Зависимость N от итерации'.decode('utf-8'))
plt.legend()
plt.grid()
plt.show()
plt.figure(figsize=(20,8))
plt.plot(np.arange(df.shape[0]), df['E'], label = 'Данные'.decode('utf-8'))
plt.xlabel('Номер итерации'.decode('utf-8'))
plt.ylabel('E'.decode('utf-8'))
plt.title('Зависимость энергии от итерации'.decode('utf-8'))
plt.legend()
plt.grid()
plt.show()
plt.figure(figsize=(20,8))
plt.plot(np.arange(df.shape[0]), df['E_r'], label = 'Данные'.decode('utf-8'))
plt.xlabel('Номер итерации'.decode('utf-8'))
plt.ylabel('E_r'.decode('utf-8'))
plt.title('Зависимость опорного потенциала от итерации'.decode('utf-8'))
plt.legend()
plt.grid()
plt.show()
E_ds = pd.DataFrame(columns=['ds', 'E'])
def analyse(frame, frameNE, ds, E_ds, begin = 500):
distr(frame)
drawNE(frameNE)
E = frameNE[begin:]['E'].mean()
s_E = np.std(frameNE[begin:]['E'])/np.sqrt(frameNE[begin:].shape[0])
E_r = frameNE[begin:]['E_r'].mean()
s_E_r = np.std(frameNE[begin:]['E_r'])/np.sqrt(frameNE[begin:].shape[0])
print E, s_E
print E_r, s_E_r
s = pd.Series([ds,E],index=['ds','E'])
E_ds = E_ds.append(s,ignore_index=True)
return E_ds
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.11/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.11/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.11, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.12/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.12/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.12, E_ds = E_ds)
del frame
del frameNE
E_ds
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.15/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.15/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.15, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.095/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.095/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.095, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.096/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.096/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.096, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.097/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.097/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.097, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.098/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.098/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.098, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.099/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.099/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.099, E_ds = E_ds)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.08/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.08/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.08, E_ds = E_ds, begin = 800)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.07/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.07/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.07, E_ds = E_ds, begin = 1000)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.06/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1000)_ds=0.06/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.06, E_ds = E_ds, begin = 1500)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x2000)_ds=0.05/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x2000)_ds=0.05/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.05, E_ds = E_ds, begin = 1500)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.04/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.04/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.04, E_ds = E_ds, begin = 3000)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.03/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x4000)_ds=0.03/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.03, E_ds = E_ds, begin = 3500)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.065/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.065/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.065, E_ds = E_ds, begin = 1000)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.075/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.075/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.075, E_ds = E_ds, begin = 1000)
del frame
del frameNE
frame = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.085/output.csv')
frameNE = pd.read_csv('D:/Backup(work)/data_output/output(1000000x1500)_ds=0.085/outputNE.csv')
E_ds = analyse(frame, frameNE, ds = 0.085, E_ds = E_ds, begin = 1000)
del frame
del frameNE
E_ds = E_ds.sort_values(by=['ds'])
E_ds.to_csv('E_ds.csv')
plt.plot(E_ds['ds'], E_ds['E'])
plt.grid()
plt.show()
E_ds